Stochastic Acceleration in Relativistic Parallel Shocks 



O 
O 

(N 

> 
O 

00 



Joni J. P. Virtanen 

Tuorla Observatory 1 
Vdisdldntie 20, FI-21500 Piikkid, Finland 
j oni . virtanenOutu . f i 
and 

Rami Vainio 

Department of Physical Sciences 
P.O. Box 64, FI-00014 University of Helsinki, Finland 
rami . vainio@helsin.ki . f i 



> 
OO 



o 



Or 

6 

CO 

C3 



X 



ABSTRACT 

We present results of test-particle simulations on both the first and the second order Fermi 
acceleration (i.e., stochastic acceleration) at relativistic parallel shock waves. We consider two 
scenarios for particle injection: (i) particles injected at the shock front, then accelerated at the 
shock by the first order mechanism and subsequently by the stochastic process in the downstream 
region; and (ii) particles injected uniformly throughout the downstream region to the stochastic 
process mimicing injection from the thermal pool by cascading turbulence. We show that re- 
gardless of the injection scenario, depending on the magnetic field strength, plasma composition, 
and the employed turbulence model, the stochastic mechanism can have considerable effects on 
the particle spectrum on temporal and spatial scales too short to be resolved in extragalactic 
jets. Stochastic acceleration is shown to be able to produce spectra that are significantly flatter 
than the limiting case of N(E) oc E^ 1 of the first order mechanism. Our study also reveals a 
possibility of re-acceleration of the stochastically accelerated spectrum at the shock, as particles 
at high energies become more and more mobile as their mean free path increases with energy. 
Our findings suggest that the role of the second order mechanism in the turbulent downstream 
of a relativistic shock with respect to the first order mechanism at the shock front has been 
underestimated in the past, and that the second order mechanism may have significant effects on 
the form of the particle spectra and its evolution. 

Subject headings: acceleration of particles — shock waves — turbulence 



1. INTRODUCTION 

For particle acceleration in astrophysical sources, 
such as jets and shock waves, two kinds of Fermi 
acceleration mechanisms have typically been con- 
sidered: the first order Fermi acceleration at the 
shock fronts, and the second order Fermi acceler- 
ation (often referred to as stochastic acceleration) 
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in the turbulent plasma. The first order mecha- 
nism is well known to produce a power-law particle 
spectrum N(E) oc E~ s with spectral index being 
a function s = (r + 2)/(r — 1) of the compres- 
sion ratio r for non-relativistic shocks (e.g., Drury 
1983) and in relativistic shocks approaching the 
value s s» 2.2 at the ultra-relativistic limit (e.g., 
Kirk et al. 2000; Achterberg et al. 2001; Lemoine 
& Pelletier 2003; Virtanen & Vainio 2003a), and 
it has been successfully applied in explaining the 
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synchrotron properties of, e.g., some active galac- 
tic nuclei (AGN) (e.g., Tinier et al. 2000). For 
flatter spectra (s < 2) the first order mechanism 
is not, however, as attractive. Although it can 
produce hard spectra with spectral index s — ► 1 
depending, for instance, on the injection model 
or the scattering center compression ratio (e.g., 
Ellison ct al. 1990; Vainio, Virtanen & Schlick- 
eiser 2003, hereafter referred to as VVS 2003), it 
is not able to produce spectral indices flatter than 
s = 1. Stochastic acceleration, on the other hand, 
has been known to exist and be present in the 
turbulent downstream of shocks, but because it 
works on much longer time scales than the first 
order mechanism (e.g., Campeanu & Schlickciscr 
1992; Vainio & Schlickeiser 1998) it has been fre- 
quently neglected in studies of relativistic parti- 
cle acceleration (however, note e.g., Ostrowski & 
Schlickeiser 1993). 

The argument of longer acceleration time scale 
renders the second order mechanism less impor- 
tant than the first order one, when the two mecha- 
nisms operate concurrently and, thus, the particle 
spectrum at the shock front is considered. When 
discussing non-thermal particle distributions radi- 
ating in astrophysical objects, however, it is im- 
portant to acknowledge that the bulk of radiation 
is emitted by the particles that have already left 
the shock front towards the downstream. Thus, 
the second order mechanism has a much longer 
time available to accelerate the particles than the 
first order mechanism. Although one could justify 
the neglect of stochastic acceleration when calcu- 
lating the particle spectrum right at the shock 
front, it is not possible to neglect its effect on 
the spectrum of radiating particles in astrophys- 
ical shock waves in general as these objects, es- 
pecially those believed to host relativistic shock 
waves, are often spatially unresolved. 

In this paper we have studied the possible ef- 
fects of stochastic electron acceleration in parallel 
relativistic shock waves. We approach the sub- 
ject via numerical test-particle simulations and 
present our model including both the first and the 
second order Fermi acceleration. We employ the 
model to study the effect of stochastic accelera- 
tion on (i) particles injected at the shock front and 
subsequently accelerated by the first order mech- 
anism and (ii) particles drawn from the heated 
(but not shock accelerated) particle population of 



the downstream region of the shock. We focus on 
shock waves that, in addition to being parallel, 
have small-to-intermediate Alfvenic Mach num- 
bers. Low Mach-number relativistic shocks could 
prevail in magnetically dominated jets that are 
lighter than their surroundings, e.g., because of 
having a pair-plasma composition. 

The structure of this paper is as follows. In 
§2 we present our model, state the limiting as- 
sumptions we have made, and briefly discuss the 
limitations caused by the assumptions and simpli- 
fications; the description of the numerical code, as 
well as the implementations and mathematical de- 
tails of the underlying physics are described in Ap- 
pendix A. In §3 we present the results achieved us- 
ing the simulation code. We begin by showing that 
when compared to some relevant previous studies 
- both analytical and numerical - our results are 
in very good accordance with those achieved pre- 
viously by many authors. We then continue pre- 
senting the results for stochastic acceleration in 
various cases, and in §4 we discuss our results and 
their relationship to both the previous studies and 
possible future applications, and list the conclu- 
sions of our study. Also the limitations caused by 
the test-particle approach are shortly discussed. 

2. MODEL 

In this section we describe the properties of our 
model in general, including the employed assump- 
tions and physics related to them. The numeri- 
cal Monte-Carlo approach is described in detail in 
Appendix A, where also implementations of the 
model properties are discussed. 

2.1. Coordinate systems 

Before proceeding into the model, we define the 
coordinate systems employed in our study: the 
frame where the shock front is at rest is called 
the shock frame. We consider parallel shocks, 
i.e., shocks where the mean magnetic field and 
the plasma flow are directed along the shock nor- 
mal. The frame where the bulk plasma flow is at 
rest is called the local plasma frame; this frame 
moves with the local flow speed, V, with respect 
to the shock frame. Finally, we consider the wave 
frame, which moves with the phase speed of 
the plasma waves with respect to the local plasma 
frame and, thus, at speed (V + V<t,)/{1 + VV^/c 2 ) 
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with respect to the shock frame denoting, as usual, 
the speed of light by c. If the scattering centers are 
taken to be fluctuations frozen-in to the plasma 
then the speed of the waves with respect to the 
underlying plasma flow is V^, = and the plasma 
frame is also the rest frame of the scattering cen- 
ters. 

2.2. Shock Structure 

We use the hyperbolic tangent profile of Schnei- 
der & Kirk (1989) to model the flow velocities at 
different distances from the shock. The width of 
the transition from the far upstream flow speed, 
V\, to that of the far downstream, V2, takes place 
over a distance of 0.01A e (ri) (for which the shock 
can still be considered almost step- like; see, e.g., 
Virtanen & Vainio 2003a), where A e (7) denotes 
the mean free path of the electrons as a function 
of Lorcntz factor 



7=1/^1-^ 



(1) 



v is the electron speed, and Ti is the Lorentz factor 
of the upstream bulk flow. (We use the standard 
notation of subscript 1[2] denoting the upstream 
[downstream] values.) The ratio of the flow speeds 
on both sides of the shock gives the gas compres- 
sion ratio 

We fix the shock speed in the upstream rest frame 
and, thus, the (equal) upstream bulk speed V\ in 
the shock frame. Using this we compute the corre- 
sponding gas compression ratio for a given shock 
proper speed 
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(3) 



following a scheme described by VVS (2003) and 
shown in Figure 1. Finally, the flow speed in the 
far downstream V2 is given by equation (2). 

2.3. Magnetic Field and Scattering 

For modeling the magnetic field structure we 
adopt the quasilinear approach, where the field B 
is considered to consist of two parts: the static, 
large scale background field B , and smaller scale 
fluctuations with amplitude SB < B n . We model 
the turbulence as being composed of a wide-band 
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Fig. 1. — Gas compression ratio r in a parallel 
shock propagating into a cold medium with proper 
speed ui (VVS 2003). 

spectrum of Alfvcn waves propagating both par- 
allel and anti-parallel to the flow. In the local 
plasma frame the waves propagate at Alfven speed 
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where h = (p + P) / n = w/n, n, p, and P are the 
specific enthalpy, the number density, the total en- 
ergy density, and the gas pressure - all measured 
in the local plasma frame. The wave intensity as a 
function of wavenumber, I(k), is assumed to have 
a power-law form for wavenumbers above an in- 
verse correlation length fco- We write the power- 
law as 

I(k)=I {k /ky, k>k (5) 

where q is the spectral index of the waves. For 
wavenumbers smaller than fco the wave intensity 
per logarithmic bandwidth is assumed to be equal 
to the background field intensity, i.e., I(k) = 
Bq^ 1 for k < ko- In this work we use two val- 
ues for q: 2 and 5/3. The former produces rigid- 
ity independent scattering mean free paths, while 
the latter is consistent with the Kolmogorov phe- 
nomenology of turbulence. 

Electrons scatter off the magnetic fluctuations 
resonantly. The scattering frequency of electrons 
with Lorentz factor 7 is determined by the inten- 
sity of waves at the resonant wavenumber (see Ap- 
pendix A.) 
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is the relativistic electron gyrofrequency, and 
^e,o = is its non-relativistic counterpart. 

Scatterings are elastic in the wave frame and 
the existence of waves propagating in both direc- 
tions at a given position, thus, leads to stochastic 
acceleration (Schlickeiser 1989). Since the spec- 
trum of waves is harder below k = k , scattering 
at energies 



7 > 7o 
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becomes less efficient. Thus, we expect the elec- 
tron acceleration efficiency to decrease at 7 > 70. 
Instead of trying to fix the value of fco, we use a 
constant value 70 = 10 6 , which corresponds to ob- 
servations of maximum Lorentz factor of electrons 
in some AGN jets (Meisenheimer et al. 1996). 

In addition to scattering, the particles are also 
assumed to lose energy via the synchrotron emis- 
sion. The average rate of energy loss for an elec- 
tron with Lorentz factor 7 in the frame co-moving 
with the plasma is given by 



dE 4 2 
_ = --a T c 7 U B , 



(9) 



where <tt is the Thompson cross-section and 
Ub = Bq/(8it) is the magnetic field energy den- 
sity. We calculate the latter in all simulations by 
assuming a hydrogen composition of the plasma. 

2.4. Alfven-Wave Transmission 

Downstream Alfven-wave intensities can be 
calculated from know upstream paramters (e.g., 
Vainio & Schlickeiser 1998; VVS 2003). Regard- 
less of the cross helicity of the upstream wave 
field (only parallel or anti-parallel waves, or both), 
there are always both wave modes present in the 
downstream region. The transmission coefficients 
for the magnetic field intensities of equation (5) 
of forward (+) and backward (— ) waves at con- 
stant wavenumber k, measured in the wave frame, 
can be written (see VVS 2003, equations (22) and 
(23)), as 
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Using these we can solve the amplification factor 
of the wave intensity for both wave modes for con- 
stant wave frame wave number k as 



W- k± =f? ± +Rl T . 



(11) 
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Fig. 2. — Ratio of the amplified wave intensities as 
function of Alfvenic Mach number. Cases corre- 
sponding to proper Alfvcn speeds ua,i of 10.0c, 
1.0 c, and 0.1c are plotted with solid, dotted, 
and dashed lines, respectively, in the case of Kol- 
mogorov turbulence, to q = 5/3. 

Amplification factor depends on the strength of 
the magnetic field as well as on the form of the tur- 
bulence spectrum. The intensity ratio of the back- 
ward waves to that of forward waves as a function 
of the quasi-Newtonian Alfvenic Mach number, 



M = Ui/« A ,1, 



(12) 



is plotted in Figure 2 for three different Alfvcn 
proper speed 
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where (3a,i = «A,i/ c - The waves are seen to prop- 
agate predominantly backward for relatively low- 
Mach-number shocks as shown by VVS (2003) in 
case of relativistic shocks, and by Vainio & Schlick- 
eiser (1998) for non-relativistic shocks. This en- 
ables the scattering center compression ratio Tk to 
grow larger than the gas compression ratio r and, 
thus, to cause significantly harder particle spectra 
compared to the predictions of theories relying on 
fluctuations frozen-in to plasma flow. As the Mach 
number increases, the downstream wave intensi- 
ties approach equipartition at the ultra-relativistic 
limit. 

Waves conserve their shock-frame frequencies 
during the shock crossing (Vainio & Schlick- 
eiser 1998), and for given upstream wave- frame 
wavenumbers k\± (here equipartition of the up- 
stream waves is assumed) the downstream wave 
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frame values are obtained from (VVS 2003) 



T 1±Vl± 

i 2+ V2+ 



Tv^7 fcl± (14) 



for both the forward (&2+) and the (&2-) waves. 
Here V}± and Tj± refer to the wave speeds of for- 
ward (+) and backward (— ) waves in the upstream 
(j = 1) and downstream (j = 2) region and to the 
respective Lorentz factors. The functional form of 
the spectrum does not change on shock crossing 
and the spectral index q in equation 5 is the same 
both sides of the shock. 

2.5. Testing the Model 

We have tested the ability of our model to 
produce results expected from previous numer- 
ical studies and theory. To test the model in 
the case of first order acceleration we ran numer- 
ous test-particle simulations with different injec- 
tion energies and shock widths (Virtanen & Vainio 
2003a, b). We found our results to be in very 
good agreement with the semi-analytical results 
for modified shocks of Schneider & Kirk (1989), 
as well as to the numerical results of, e.g., Elli- 
son et al. (1990) for the corresponding parts of 
the studies. For the step-shock approximation, 
spectra with indices close to the predicted value of 
<~ 2.2 were obtained. An example of the test-runs 
is shown in Figure 3 with the shock proper speed 
u\ = 10 c and compression ratio r « 3, assuming 
scattering centers frozen-in into the plasma and 
turbulence having a spectral index of q = 2. 

For testing the second order acceleration we 
chose an analytically well known case, namely that 
of uniform flow with waves streaming in parallel 
and anti-parallel directions with equal intensities. 
In this case the momentum diffusion coefficient 
of charged particles can be given as (Schlickeiser 
1989) 

2^ c 2 -%% p 2 v\ 



A 2 (x,p) ~ 



q(q + 2)B% v 3 -l' 



(15) 



which in the relativistic case (v s=s c and p m e c) 
is found to depend on the momentum p as 



A 2 oc p q 



(16) 



For a constant particle injection at low momenta, 
this leads to a simple relation of the spectral index 
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Fig. 3. — Spectrum of the particles accelerated at 
the shock front in the step-shock approximation. 
Resulting spectral index for power-law - extending 
over five orders of magnitude in energy - is s ~ 2.2. 

of the volume-integrated particle energy spectrum, 
s, and the magnetic field fluctuations, q: 



q-l. 



(17) 



As a test case, we ran simulations involving only 
stochastic acceleration, and found that for val- 
ues qE [1,2] the model produces exactly those 
indices expected from the analytical calculation. 
The spectral indices obtained from the simulation 
are plotted together with the theoretical predic- 
tion in Figure 4. 

3. RESULTS 

In this section we apply our model to stochastic 
particle acceleration in the downstream region of 
a relativistic parallel shock, using a test-particle 
approach. 

First we show how the stochastic process affects 
a non-thermal particle population, already accel- 
erated at the shock via the first order mechanism. 
Then we consider particles injected throughout 
the downstream region. Finally we present an 
example of the combination of both injection 
schemes. Simulations were run separately for 
low, intermediate and high Alfvenic-Mach-number 
shocks (M = 3, M = 10, and M = 1000, respec- 
tively - see Fig. 2 for the corresponding wave 
intensity ratios), and for four cases of downstream 
turbulence; for turbulence spectral index q = 2 
and q = 5/3 with downstream wave field calcu- 
lated using wave transmission analysis described 
in §2.4, and with the downstream forward and 
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Fig. 4. — Comparison of the simulation model and 
the theory. Spectral index of the energy spectrum 
of particles accelerated stochastically for different 
spectral indices of the magnetic field fluctuations. 
Data from the simulations are plotted with crosses 
(statistical error of each plot is 0.02 or less), and 
the theoretical prediction of s = q — 1 is plotted 
as a straight line. 

backward waves being in equipartition. 1 The 
proper speed of the shock is set to u\ = 10 c 
in all simulations. 

3.1. Electrons Injected and Accelerated at 
the Shock 

We have studied the effect of stochastic accel- 
eration on particles that have been already accel- 
erated at the shock. This was done by injecting 
particles into the shock and the first order mecha- 
nism, and allowing them to continue accelerating 
via the stochastic process in the downstream re- 
gion. Injection of the particles took place in the 
downstream immediately after the shock, and par- 
ticles were given an initial energy of a few times 
the energy of the thermal upstream particles as 
seen from the downstream. This kind of injec- 
tion simulates some already-energized downstream 
particles returning into the shock, but without 
the need of processing the time consuming bulk 
of non-accelerating thermal particles. The high- 
energy part of the particle energy distribution 
which we are interested in in this study - is simi- 
lar, regardless of the injection energy. 

We found that in the case of high Alfvenic Mach 
number (M = 1000, corresponding to magnetic 

1 In the case of M = 1000 the effects were - at best - barely 
visible, as expected. For this reason these results are not 
included in this paper, but are available in an electronic 
form at http://www.astro.utu.fi/red/qshock.html. 



field Bq ~ 1.4 mG in a hydrogen plasma) the 
contribution of the stochastic process to the en- 
ergy distribution of the particles is, indeed, very 
insignificant compared to that of the first order 
acceleration at the shock; the energy spectrum 
maintains its shape and energy range unaltered at 
least for tens of thousands thermal electron mean 
free paths. This is the case regardless of the ap- 
plied turbulence spectral index, and because the 
analytically calculated wave intensities are very 
close to equipartition for high-Mach-number shock 
(see Fig. 2), the difference between the analyt- 
ically calculated downstream wave field and the 
explicitely assumed equipartition case was, as ex- 
pected, minimal. 

For stronger magnetic fields (M = 10 and 
M = 3, corresponding to 0.14 G and 0.46 G, re- 
spectively, in a hydrogen plasma, and to 4.6 mG 
and 15 mG in a pair plasma) the effect of stochas- 
tic acceleration is, on the contrary to the high- 
M case, very pronounced. The stochastic process 
begins to re-accelerate particles immediately after 
the shock front, and the whole spectrum slowly 
shifts to higher energies as shown in Figure 5, 
where the results obtained using wave transmis- 
sion analysis are presented in the left-hand col- 
umn, while the right-hand panels show the results 
of the equipartition assumption. The acceleration 
rate depends on the wave spectrum in a way that 
for q = 2, for which particles of all energies have 
the same mean free path, the stochastic process 
accelerates particles to higher energies at constant 
rate, whereas for the Kolmogorov turbulence the 
acceleration rate (like the mean free path of the 
particles) decreases as the energy increases, and 
the energization gradually slows down as the bulk 
of the particle energy distribution rises to higher 
energies. Also different composition of the down- 
stream wave field leads to slightly different results: 
in the "transmitted" case (left-hand panels) the 
particle population immediately behind the shock 
front extends to slightly higher energies than in the 
equipartition case (right-hand panels), but for the 
latter the stochastic acceleration is clearly quicker. 
This can be seen best in the uppermost panels of 
Figure 5, where the calculated transmitted wave 
field differed from the equipartition assumption 
the most. 

The "turnover" of the scattering rate at en- 
ergy 7o = 10 6 causes the rate of energization to 
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Fig. 5. — Evolution of the particle spectra in the turbulent downstream region due to stochastic acceleration. 
Particles are initially injected and accelerated at the shock front, located at x — near the left edge of the 
plots. Contours of log(-E^) show the steady-state particle energy distribution. On the left-hand panels the 
wave transmission calculation of VVS (2003) is applied, whereas on the right-hand panels the downstream 
waves modes are assumed to be in equipartition. In the two uppermost rows q = 2, and in the two lowermost 
panels q = 5/3. Results for Alfvenic Mach numbers M = 3 (first and third row) and 10 (second and fourth 
row) are shown for each case. The physical sizes of the downstream region are approximately 10 12 cm (first 
row), and 10 13 cm (second row) for the q — 2 cases, and 10 11 cm (third row), and 10 12 cm (fourth row) for 
the q = 5/3 cases. 
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Fig. 6. — Energy distribution of the particles es- 
caping the shock region and getting absorbed in 
the downstream. The Alfvenic Mach number of 
the shock is M = 10, the proper speed is u\ — 10 c, 
and the turbulence corresponds to q = 5/3. Note 
that the distribution is plotted as \og(dN / dE) . 

go down for energies greater than 70 . This is be- 
cause the energy dependence of the mean free path 
of the particle changes when particles resonant 
wavenumber (eqn. (6)) decreases below fc , which 
was use to set the lower limit for the I(k) oc k~ q 
power-law. For particles with 7 > 70 the mean 
free path starts to increase much faster, leading 
not only to the decrease of the stochastic accel- 
eration efficiency, but also to another notable ef- 
fect: as its scattering mean free path increases at 
7 > 70 , the particle will suddenly be able to move 
much more easily in the downstream region and 
even be able to return back to the shock. Parti- 
cles already energized first in the shock by the first 
order mechanism and then in the downstream by 
the second order acceleration, and returning back 
into the shock again, get "re-injected" into the first 
order acceleration process. This effect, in general, 
is seen in all simulations with M — 3 as a bend- 
ing of the countours to the left at 7 > 70 close 
to the shock (at x = 0), but it is visible also in 
M = 10 shocks in spectra of particles collected at 
the downstream free-escape boundary. An exam- 
ple of the latter case is shown in Figure 6. 

3.2. Electrons Injected from the Down- 
stream Bulk Plasma 

Our second approach was to assume that a con- 
stant injection mechanism exists throughout the 
downstream region and investigate the stochas- 



tic acceleration process. (Physically, this mim- 
ics a case, where turbulent fluctuations cascade 
to higher wavenumbers and inject a fraction of 
the thermal electrons to the stochastic accelera- 
tion process.) We injected particles at constant 
energy - corresponding to the energy equal to 
the energy of upstream electrons as seen from the 
downstream region - uniformly and isotropically 
within the whole downstream region. The results 
for different cases are shown in Figure 7. 

The Figure shows similar behavior as in the case 
of particles injected at the shock: significant ac- 
celeration was seen only for shocks with strong 
magnetic field, while in the M — 1000 case prac- 
tically no acceleration is seen. For the cases with 
strong magnetic field the acceleration is seen to 
work similarly as in the previous shock-related in- 
jection case, energizing particles at constant rate 
for q = 2 turbulence, and with decreasing rate for 
q = 5/3. Again the acceleration rate slows down 
when energies corresponding to 70 are reached. 
After this energy particles start to pile up and 
form a bump in the distribution immediately be- 
hind the 70. Also here (at least for M = 3) some 
particles with energy 7 > 70 are able to return 
to the shock and get re-injected to the first order 
process (see, e.g., the right-hand panel of Fig 8). 

In the energy range between the injection en- 
ergy and 70, particles begin to form a power-law 
distribution with spectral index depending on the 
magnetic field fluctuations, as expected from equa- 
tion (17). For q = 2, the particle spectral index 
s ~ 1, and for q = 5/3, s ~ 0.6. The formation 
of the power laws as function of distance from the 
shock is shown in the left-hand panel of Figure 8 
for q = 2, and in the right-hand panel for q = 5/3. 
In the latter case also the formation of high-energy 
bump at the shock by the returning accelerated 
particles is seen. 

3.3. Example of Combination 

We have investigated what kind of particle en- 
ergy spectra the two discussed injection schemes 

- one operating at the shock, and another operat- 
ing uniformly throughout the downstream region 

- are able to create. Next we will present an ex- 
ample of a combination of these. In the simu- 
lation these two cases were kept separate for the 
sake of simplicity, but there should be no reason to 
assume the separation be present also in nature. 
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Fig. 7. — Same as Fig. 5 but for particles injected uniformly and isotropically across the downstream region. 
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Fig. 8. — Energy spectrum of the particles at different distances from the shock. The particles are injected 
in the downstream region and the distributions are from the simulation with q = 2 and M = 10 (wave 
field transmitted) for the left-hand panel and from q = 5/3 and M = 3 (wave field in equipartition) for the 
right-hand panel (corresponding to the left-hand panel of the second row and the top right panel of Fig. 7). 
Slices are from locations x = (solid line), x = 100 (dashed), x = 400 (dotted), x — 1000 (dash-dotted), 
and x = 2000 (dash-dot-dot-dotted) for the left-hand panel, and from x = (solid line), x = 300 (dashed), 
x = 600 (dotted), x — 1500 (dash-dotted), and x — 3000 (dash-dot-dot-dotted) for the right-hand panel. 



Also the relative amounts of shock-injected and 
downstream-injected particles is not fixed here, 
but instead considered more or less free a param- 
eter. 

Assuming different ratios of injected particle 
populations, the resulting spectrum would be very 
different. An example of combination of the two 
injection schemes in the case of shock with u\ = 
10 c and M = 10, and with turbulence correspond- 
ing to q = 5/3, is shown in Figure 9. 

4. DISCUSSION AND CONCLUSIONS 

We have studied stochastic particle acceleration 
in the downstream region of a relativistic parallel 
shock. Applying the wave transmission calcula- 
tions of VVS (2003) and assuming the cross helic- 
ity to vanish in the upstream, we have modeled the 
turbulence of the downstream region as a super- 
position of Alfven waves propagating parallel and 
anti-parallel to the plasma flow. Using a kinetic 
Monte-Carlo simulation we have modeled the sec- 
ond order Fermi acceleration of electrons in the 
shock environment, and considered cases of accel- 
eration of downstream-injected particles, as well 
as that of particles injected at the shock. We have 
shown that the stochastic acceleration can, indeed, 
have remarkable effects in both cases. This result 



is even more pronounced if the two downstream 
Alfven wave fields are assumed to be in equiparti- 
tion. 

The behavior of the particle energy distribution 
in the stochastic process depends heavily on the 
strength of the background magnetic field; in the 
cases of weak magnetic field and quasi-Newtonian 
Alfvenic Mach number much larger than the crit- 
ical Mach number (M 3> M c = ^Jr) the effects of 
stochastic acceleration are faded out by the much 
stronger first order acceleration. Also the mag- 
netic field turbulence spectrum affects the acceler- 
ation efficiency: for Kolmogorov turbulence with 
g = 5/3 the spatial scales are up to an order of 
magnitude shorter than in the case of q = 2 turbu- 
lence. Although the spatial scales in simulations 
presented here are enormous compared to those 
associated with shock acceleration (the first or- 
der process in the immediate vicinity of the shock 
front), in case of blazars and other AGNi the scales 
are still orders of magnitude too small to be re- 
solved even in the VLB1 pictures - regardless of 
the turbulence and used magnetic field strength. 
Also the acceleration time scales are very short: 
the time required to shift the whole spectrum from 
the initial energy range to 7buik ^ 10 6 ranged from 
10 to 50 minutes in the M = 10 case, and for 
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Fig. 9. — Example of a combination (solid line) 
energy spectra of particles injected to the accelera- 
tion process at the shock (dashed) and throughout 
the downstream region (dotted). (Model parame- 
ters are q — 5/3 and M = 10.) The particles are 
collected at the downstream free escape boundary, 
~4x 10 12 cm away from the shock, and the num- 
ber of the particles injected at the shock is larger 
than the number of uniformly injected particles by 
a factor of 100. 

M = 3 the times were < 1 minute, as measured in 
the shock frame. 

In addition to the magnetic field strength and 
turbulence, also the composition of the down- 
stream wave field was seen to affect the result- 
ing particle population. When comparing other- 
wise similar cases that differ only for the down- 
stream cross helicity (i.e., whether the wave field 
is resulting from the wave transmission calcula- 
tions of VVS (2003) or an equipartition of parallel 
and anti-parallel waves is assumed), the calculated 
wave-transmission cases with more anti-parallel 
waves (VVS 2003, and Fig. 2) showed stronger first 
order acceleration, but weaker stochastic accel- 
eration. This is because of the larger scattering 
center compression ratio in the wave-transmission 
case leading to more efficient first order acceler- 
ation (VVS 2003) and, on the other hand, faster 
momentum diffusion rate in the equipartition case 
leading to more efficient stochastic acceleration. 

In the cases where the stochastic acceleration 
was quick enough for particles to reach the 70 en- 
ergy while being still sufficiently close to the shock 
in order to be able to make their way back to the 



upstream region due to their prolonged mean free 
path, the first order mechanism was able to re- 
accelerate the returning high-energy particles to 
even higher energies. This led to forming of a 
new (quasi-)power-law at energies 7 > 10 7 in some 
cases. 

One notable feature of the present model is that 
in the case of a uniform injection process in the 
downstream region, power-law spectra with high 
and low energy cut-offs are formed. Depending 
on the turbulence, particle energy spectra have 
power-law spectral indices of 0.5-1 with lower and 
higher energy cut-offs at 71 w 10 1 ~ injection and 
72 ~ 10 6 = 70, respectively. These particles would 
produce synchrotron spectra with photon spectral 
indices — 0.5<a<0in the GHz-THz regime for 
various initial parameters. These properties are 
quite similar to those of flat-spectrum sources, for 
which typical spectra with a > —0.5 in the GHz 
region and flare spectra with a w —0.2 in the op- 
tically thin region of the spectrum are seen (e.g., 
Valtaoja et al. 1988, and references therein). 

Combining the resulting energy distributions of 
both the particles injected at the shock and the 
particles injected uniformly in the downstream re- 
gion can lead to very different results depending on 
the relative amounts of particles injected in both 
cases. Because different distributions produce var- 
ious observable spectra, one might, in principle, be 
able to set constraints on different injection mech- 
anisms, as well as the physical size of the tur- 
bulent downstream region of the shock, compar- 
ing the evolution of observed radiated spectra and 
the predictions basing on several composite distri- 
butions. Especially the maximum Lorentz factor 
of electrons could be used to set limits for 70, as 
well as the lowest-energy parts of the non-thermal 
power-law spectra could give hints of the injection 
process. 

Our simulations are based on test-particle ap- 
proximation, i.e., the effects of the particles on the 
turbulent wave spectrum and on the shock struc- 
ture are neglected. Including wave-particle inter- 
actions in a self-consistent manner may modify the 
cross-helicities and wave intensities in the down- 
stream region and lead to notable effects on the 
accelerated particle spectrum (e.g., Vainio 2001). 
One should note, however, that the wave particle 
interactions are competing with wave-wave inter- 
actions in the turbulent downstream region, which 
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may modify the the turbulence parameters in a 
different manner. Including these effects to our 
model are, however, beyond of the scope of the 
present simulations. 

To conclude, the main results of this paper are: 
(i) Stochastic acceleration can be a very efficient 
mechanism in the downstream region of paral- 
lel relativistic shocks, provided that the magnetic 
field strength is large enough in order to make the 
Alfvenic Mach number approach the critical Mach 
number (M c = y/r) of the shock, i.e., to increase 
the downstream Alfven speed enough to allow for 
sufficient difference in speeds of parallel and anti- 
parallel Alfven waves required for rapid stochastic 
acceleration, (ii) The forming of a power-law with 
very hard particle energy spectra between the in- 
jection energy and 70 in the case of a continuous 
injection mechanism in the downstream region. 
The produced particle populations could produce 
synchrotron spectra very similar to those of flat 
spectrum sources, (iii) The interplay between the 
first and second order Fermi acceleration at rel- 
ativistic shocks can produce a variety of spectral 
forms not limited to single power laws. 

J.J.P.V. thanks the Finnish Cultural Foun- 
dation for financial support. 

Appendix A. The Monte Carlo Code 

In this appendix we review the structure and 
implementation of our simulation code. The code 
employs a kinetic test-particle approach; it follows 
individual particles in a pre-defined and simplified 
shock environment, based on the assumptions and 
simplifications presented in § 2. In short, the sim- 
ulation works as follows: we trace test particles 
under the guiding center approximation in a ho- 
mogeneous background magnetic field with super- 
posed magnetic fluctuations. The fluctuations are 
assumed to be either static disturbances frozen-in 
into the plasma flow, or Alfven waves propagating 
along the large-scale magnetic field parallel and 
antiparallel to the flow (in this paper we apply 
only the Alfven wave case). In each time-step the 
particle is moved and scattered; scatterings are 
modeled as pitch-angle diffusion with an isotropic 
diffusion coefficient cx 1 — /i 2 , where \i = cos 9 
is the pitch-angle cosine of the particle. Also for 
each simulation time-step the energy losses due to 



the synchrotron emission (see cq. (9)) are calcu- 
lated. 

When the particle passes a free-escape bound- 
ary in the far downstream region it is removed 
from the simulation. The injection of the parti- 
cles into the simulation is modeled using several 
methods. The first method involves a uniform and 
isotropic injection of particles within one mean 
free path of a thermal electron downstream of 
the shock front, simulating the already-energized 
supra-thermal particles crossing and re-crossing 
the shock. This injection method allows us to 
concentrate on the non-thermal particles without 
having to spend most of the computing time sim- 
ulating the thermal body of the total particle dis- 
tribution. Other injection employed methods in- 
clude an injection of a cold distribution upstream, 
or an injection of particles at thermal energies uni- 
formly and isotropically in the downstream region 
(the latter case is applied in §3.2). 

The time unit of the simulation is chosen to 
be the inverse of the scattering frequency of an 
electron having energy Er 1 = Tim c c 2 , where Ti 
is the Lorentz factor of the shock. The unit of 
velocity is chosen as c. With these choices the unit 
of length equals the mean free path of electron 
with Lorentz factor that of the shock: A e (Fi) = 
l/v c (Ti). The time-step is chosen so that it is a 
small fraction of the inverse scattering time, At = 
av c {"f), where a < 0.01. 

For a typical simulation ~ 10 5 particles are in- 
jected. The number of high-energy particles is fur- 
ther increased by applying a splitting technique; if 
the energy of a particle exceeds some pre-defined 
value, the particle is replaced by two "daughter 
particles" which are otherwise identical to their 
"mother" , but have their statistical weight halved. 
The number of these splitting boundaries is chosen 
so that for each simulation the balance between 
the statistics and the simulation time is optimal. 

A.l. The Shock and the Flow Profile 

We consider a shock wave propagating at speed 
V\ into a cold ambient medium, which is taken to 
be initially at rest in the observer's frame. In the 
shock frame the shock is, by definition, at rest and 
the upstream plasma flows in with speed V\. At 
the shock front the plasma slows down, undergoes 
compression, and flows out at downstream flow 
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speed V 2 < V\. The shocked downstream values 
of the plasma parameters are determined by the 
compression ratio of the shock defined in equation 
(2). 

For the flow speed depending on the location x, 
we use the hyperbolic tangent form of Schneider 
k Kirk (1989), 



V(x) = V 1 - 



V 1 -V 2 



1 + tanh 



WAe(ri) 



(18) 

and fix the shock thickness to W = 0.01, cor- 
responding to a nearly step-like shock. For 
thicker shocks the first order acceleration efficiency 
drops fast (e.g. Schneider & Kirk 1989; Virta- 
nen & Vainio 2003b), and also the simple wave- 
transmission calculations of VVS (2003) would 
not be valid. 

A. 2. Magnetic Field 

The strength of the parallel magnetic field, Bo, 
in the simulation is defined with respect to the 
"critical strength" for which the Alfven speed in 
the downstream region becomes equal to the local 
flow speed; for fields stronger than this the paral- 
lel shock becomes non-evolutionary. The critical 
field, B c , is calculated from equation 



B c c/^Airh 2 n 2 + B 2 c = V 2 , 



(19) 



for which we get the downstream particle density 
n 2 and the specific enthalpy h 2 using the magne- 
tohydrodynamical jump conditions (e.g., Kirk & 
Duffy 1999) and the upstream values m = 1 cm~ 3 
and hi = {pi+Pi)/n\ = mc 2 with m = m c +m p in 
a hydrogen plasma and m — 2m in a pair plasma 
(throughout this work the upstream plasma is as- 
sumed to be cold so that its pressure may be ne- 
glected compared to its rest energy) . Applying the 
equations for conservation of both energy 



h l T 1 V l = h 2 T 2 V 2 



(20) 



and mass 



riVim = T 2 V 2 n 2 , 



(21) 

and a bit of straightforward algebra, the critical 
field can be written as 



B c = ^A-kVi V 2 mni IY 



(22) 



For example, for a shock with compression ratio 
r = 3 propagating into a cold ambient medium 



of ni = 1 cm~ 3 at proper speed u\ = 10 c the 
critical field is B c ss 0.8 G for a hydrogen plasma 
and 0.03 G for a pair plasma. 

A. 3. Particle Transport and Scattering 

During each Monte Carlo time-step scatterings 
off the magnetic fluctuations are simulated by 
making small random displacements of the tip of 
the particle's momentum vector, and the parti- 
cle is transported according to its parallel (to the 
flow) speed in the fixed shock frame. In the case 
of stochastic acceleration, the particle is scattered 
twice for each Monte Carlo time-step. Instead of 
only one scattering off scattering centers frozen-in 
into the plasma, the scattering process is now di- 
vided into two parts: the particle is first scattered 
in the forward-wave frame, and then immediately 
again in the backward- wave frame. The scattering 
frequencies of both scatterings are adjusted so that 
the total statistical effect of the duplex scattering 
is consistent with the value of the momentary scat- 
tering mean free path. This is done by balancing 
the both scattering frequencies by weight values 
u>± € [0, 1], which are calculated from the ratio of 
the amplified waves (eq. 11 and Fig. 2), so that 
they satisfy the relation ui + + w_ = 1 . Neglecting, 
for simplicity of the simulations, the dependence 
of the scattering frequency on the particle's propa- 
gation direction, and assuming that the scattering 
is elastic in the rest frame of the scattering cen- 
ters, we can use the quasilinear theory: the scat- 
tering frequency of an electron as a function of its 
Lorentz factor (in the rest frame of the scattering 
centers, denoted here by prime) 7' is 



nCtok'I'ik') (y ; 



1 



(9-l)/2 



2 i B 2 



r 



(23) 



This equation is applicable to particles scatter- 
ing off waves with wavenumber larger than fco; for 
particles with Lorentz factor 7' > 7q the scatter- 
ing frequency decreases as v oc 7 _1 , as at these 
wavenumbers k'I(k') = B 2 is assumed. Scattering 
frequency as a function of particle's Lorentz factor 
is plotted in Figure 10 for magnetic field fluctua- 
tions with q = 2 (leading to energy-independent 
mean free path) and q = 5/3. 

In each scattering, the velocity vector of the 
particle is first Lorentz-transformed into the frame 
of the scatterers, then the new pitch-angle cosine 
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Fig. 10. — Scattering frequency as a function of 
Lorentz factor. The turn-over energy (eq. (8), see 
text for details) 70 = 10 6 is marked with a vertical 
dotted line. The values of u{fj) corresponding to 
q = 5/3 and q = 2 are plotted with solid and 
dashed curves, respectively. 
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Fig. 11. — The geometry of the scattering event. 
Particle is initially moving at velocity v ; the angle 
between the velocity and magnetic field B is 9q. 
(Thus, 80 = arccos^Q.) Scattering changes the 
direction of the particle by an angle i9 so that the 
new velocity v has an angle 9 between it and the 
magnetic field. <p is the azimuth angle of the new 
velocity measured from the plane defined by B and 
v - 

(in the scatterer frame) is computed from the for- 
mula (e.g., Ellison ct al. 1990) 

(j, 1 <— p! cos ■& + y/l - ii! 2 sin -Q cos 0, (24) 

where the angle between the velocities before (v ) 
and after (v) the scattering, # <G [0,7r), and the 
angle measured around the scattering axis, 4> € 
[0, 2ir) (angles and the geometry is sketched in Fig. 
11), arc picked via a random generator from ex- 
ponential and uniform distributions, respectively 
(see Vainio ct al. 2000, for details). The new veloc- 
ity vector is then Lorentz transformed back to the 
shock frame, and finally the the particle is moved 
according to its new parallel velocity. 
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